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ABSTRACT 

A  variational  data  assimilation  algorithm  is  developed  for  the  ocean  wave  prediction  model  [Wave  Model 
(WAM)].  The  algorithm  employs  the  adjoint-free  technique  and  was  tested  in  a  series  of  data  assimilation 
experiments  with  synthetic  observations  in  the  Chukchi  Sea  region  from  various  platforms.  The  types  of 
considered  observations  are  directional  spectra  estimated  from  point  measurements  by  stationary  buoys, 
significant  wave  height  (SWH)  observations  by  coastal  high-frequency  radars  (HFRs)  within  a  geographic 
sector,  and  SWH  from  satellite  altimeter  along  a  geographic  track.  Numerical  experiments  demonstrate 
computational  feasibility  and  robustness  of  the  adjoint-free  variational  algorithm  with  the  regional  configu¬ 
ration  of  WAM.  The  largest  improvement  of  the  model  forecast  skill  is  provided  by  assimilating  HFR  data 
(the  most  numerous  among  the  considered  types).  Assimilating  observations  of  the  wave  spectrum  from  a 
moored  platform  provides  only  moderate  improvement  of  the  skill,  which  disappears  after  3  h  of  running 
WAM  in  the  forecast  mode,  whereas  skill  improvement  provided  by  HFRs  is  shown  to  persist  up  to  9  h.  Space- 
borne  observations,  being  the  least  numerous,  do  not  have  a  significant  impact  on  the  forecast  skill  but  appear 
to  have  a  noticeable  effect  when  assimilated  in  combination  with  other  types  of  data.  In  particular,  when 
spectral  data  from  a  single  mooring  are  used,  the  satellite  data  are  found  to  be  the  most  beneficial  as  a 
supplemental  data  type,  suggesting  the  importance  of  spatial  coverage  of  the  domain  by  observations. 


1.  Introduction 

In  the  last  decades,  significant  progress  has  been  made 
in  numerical  modeling  of  oceanic  waves  (Cavaleri  et  al. 
2007;  Janssen  2008).  With  the  advent  of  massive  obser¬ 
vations  of  sea  surface  roughness  from  satellites,  the 
forecast  skill  of  the  wave  models  (WMs)  acquired  new 
prospects  for  further  improvement  through  optimiza¬ 
tion  of  the  poorly  known  model  parameters  with  respect 
to  the  observations. 

The  early  attempts  to  improve  WM  performance 
through  data  assimilation  employed  the  nudging  tech¬ 
nique  to  reduce  the  misfit  between  the  integral  charac¬ 
teristics  of  the  model  spectra,  such  as  significant  wave 
height  (SWH),  and  observations  (Francis  and  Stratton 
1990;  Bauer  et  al.  1992;  Breivik  and  Reistad  1994; 
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Wittmann  and  Cummings  2004).  More  sophisticated  as¬ 
similation  schemes  updated  the  model  using  optimal  in¬ 
terpolation  of  several  spectral  characteristics  (Lionello 
et  al.  1992;  Dunlap  et  al.  1998;  Greenslade  2001).  The 
approach  was  generalized  to  the  spectral  partitioning 
technique  (Voorrips  et  al.  1997;  Aouf  and  Lefevre  2006; 
Abdalla  et  al.  2006),  which  appeared  to  provide  a  rea¬ 
sonably  good  improvement  of  the  forecast  skill. 

An  interesting  approach  was  proposed  by  Bauer  et  al. 
(1996),  who  explicitly  inverted  the  linearized  Wave 
Model  (WAM;  WAMDI  Group  1988;  Monbaliu  et  al. 
2000)  operator  to  obtain  corrections  to  the  wind  fields 
induced  by  observations  of  the  wave  spectra.  This 
method  was  constrained  by  a  rather  strong  assumption 
that  the  response  of  the  linearized  model  is  localized 
both  in  space  and  time.  Holthuijsen  et  al.  (1997) 
explored  a  similar  technique  using  a  limited  number  of 
nonlocal  perturbations  that  were  derived  from  statistical 
analysis  of  the  wind  forcing  error  fields.  The  method 
demonstrated  a  20%-25%  reduction  of  the  forecast 
errors  at  time  scales  less  than  a  day. 
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Voorrips  et  al.  (1999)  were  among  the  first  to  explore 
the  performance  of  the  extended  Kalman  filter  (KF)  in 
application  to  WAM  using  simulations  in  simple  ge¬ 
ometries.  The  extended  KF,  its  reduced-order  ap¬ 
proximation,  and  fixed-lag  Kalman  smoother  were 
tested  and  they  demonstrated  better  performance  than 
optimal  interpolation  (Hasselmann  et  al.  1997),  espe¬ 
cially  in  estimating  the  wave  field  away  from  observa¬ 
tions  and  in  past  wind  corrections.  More  recently, 
Sannasiraj  et  al.  (2006)  proposed  a  hybrid  assimilation 
scheme  that  combined  a  statistical  error  model  at  the 
observation  points  with  spatial  interpolation  using  the 
Kalman  gain  matrix.  Their  method  was  capable  of  re¬ 
ducing  the  forecast  errors  by  30%-60%  over  12-24-h 
forecast  periods.  Operationally  oriented  sequential  KF 
methods  were  also  considered  by  Emmanouil  et  al. 
(2010, 2012),  who  employed  WAM  in  combination  with 
the  second-order  KF  and  optimal  interpolation 
schemes  to  improve  the  model’s  performance  in  the 
North  Atlantic. 

Variational  methods  involving  the  adjoint  models 
have  been  in  use  in  WM  data  assimilation  for  more 
than  two  decades  (Snyder  et  al.  1992;  Barzel  and  Long 
1994;  De  Las  Heras  et  al.  1994,  1995;  Hersbach  1997; 
Veeramony  et  al.  2010).  One  of  the  difficulties  with 
this  approach  is  in  the  numerical  complexity  of  the 
WMs,  which  employ  sophisticated  nonlocal  parame- 
terizations  of  the  source  terms  and  operate  in  spatially 
variable  spectral  space.  De  Las  Heras  et  al.  (1994)  used 
an  approximation  to  the  adjoint  WAM,  neglecting 
spatial  derivatives  in  the  evolution  equation.  Orzech 
et  al.  (2013)  took  an  alternative  approach  of  de¬ 
veloping  the  adjoint  for  spectral  advection  and 
performed  a  series  of  experiments  controlling  steady- 
state  open  boundary  conditions  of  the  linearized  ver¬ 
sion  of  spectral  WMs  for  nearshore  [Simulating  Waves 
Nearshore  (SWAN)]  applications.  Similar  experi¬ 
ments  were  performed  by  Veeramony  et  al.  (2010), 
who  used  the  analytical  adjoint  developed  by  Walker 
(2006).  Recently,  Orzech  et  al.  (2014)  employed  the 
adjoint  of  SWAN  to  optimize  the  offshore  SWH 
observations. 

In  the  present  study  we  explore  the  feasibility  of  the 
new  variational  data  assimilation  technique  (Yaremchuk 
et  al.  2009;  Trevisan  et  al.  2010)  in  application  to  WAM. 
This  approach  employs  a  sequence  of  low-dimensional 
subspaces  that  are  iteratively  updated  in  the  process  of 
finding  a  cost  function  minimum.  A  distinctive  feature 
of  the  method  is  that  it  does  not  require  development 
of  the  tangent  linear  and  adjoint  codes  for  imple¬ 
mentation,  which  is  important  for  WM  applications.  In 
addition,  the  adjoint-free  four-dimensional  variational 
data  assimilation  (a4DVAR)  method  is  less  vulnerable 


Fig.  1.  Wind  speed  (white  arrows)  and  SWH  (contours,  m)  of  the 
model  solution  at  0000  UTC  20  Sep  2011  ( t  =  0).  Mooring  positions 
are  shown  by  black  squares.  HFR  location  and  coverage  area  are 
given  by  the  black  circle  with  a  sector.  SWH  data  are  acquired 
along  the  radar  beams  shown  by  dotted  lines  within  the  sector. 
Dashed  lines  are  the  tested  tracks  of  Enviscit. 


to  the  instabilities  that  develop  in  the  adjoint  versions 
of  the  numerical  models  due  to  the  breakdown  of  the 
tangent  linear  approximation  in  strongly  nonlinear  re¬ 
gimes.  It  was  also  shown  that  the  a4DVAR  technique 
appears  to  be  advantageous  when  observations  are 
sparse  and  noisy. 

We  explore  the  performance  of  the  a4DVAR  tech¬ 
nique  with  WAM  configured  in  the  800  km  X  800  km 
region  of  the  Chukchi  Sea  (Fig.  1).  In  recent  years  this 
sector  of  the  Arctic  Ocean  was  subject  to  significant 
summer/fall  storm  activity  due  to  the  rapid  loss  of  the  ice 
cover,  but  its  wave  conditions  remain  largely  un¬ 
explored.  The  choice  of  the  region  was  partly  motivated 
by  the  anticipated  increase  of  shipping  through  the  Be¬ 
ring  Strait.  This  may  require  monitoring  and  forecasting 
of  the  wave  conditions,  which  could  be  better  done  with 
data  assimilation.  Presented  data  assimilation  (DA) 
experiments  are  performed  with  several  types  of  ob¬ 
servations  that  include  SWH  data  from  satellites,  coastal 
HF  radars,  and  spectral  characteristics  obtained  from 
in  situ  observations  at  moored  buoys. 

In  the  following  section,  we  present  the  details  of  ex¬ 
perimental  setting  after  a  brief  description  of  the 
a4DVAR  assimilation  technique.  In  section  3,  the  re¬ 
sults  of  synthetic  data  experiments  are  presented  and 
compared  with  the  results  of  the  traditional  optimal  in¬ 
terpolation  (OI)  method,  observability  of  the  wave  field 
from  different  platforms  is  discussed,  and  the  compu¬ 
tational  cost  of  the  a4DVAR  is  considered.  Additional 
discussion  is  offered  in  section  4,  preceded  by  a  brief 
summary  and  overall  conclusions. 
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2.  Methodology 


where 


a.  a4DVAR  technique 

A  variational  DA  algorithm  performs  a  dynamically 
constrained  search  for  a  minimum  of  the  cost  function  /, 
which  quantifies  the  distance  between  observations  and 
the  numerical  model  solution  used  to  assimilate  the  data. 
This  minimization  is  performed  with  respect  to  a  set  of 
poorly  known  model  parameters  (e.g.,  the  set  of  gridpoint 
values  of  the  model  fields  at  the  start  of  time  integration) 
that  constitute  the  control  vector  c.  Assuming  the  dif¬ 
ferentiability  of  J  with  respect  to  c,  in  the  vicinity  of  the 
minimum  the  cost  function  is  represented  by 


-B-t/2- 

■  0  ■ 

s  = 

HjM1 

;  b  = 

dl 

H  M” 

L  n  -1 

A- 

(6) 


and  S  is  the  square  root  of  the  Hessian  matrix  J  =  STS. 
To  find  the  minimum,  one  has  to  solve  the  normal 
equation 


STSc  =  STb. 


(7) 


/  =  /o+t(c~c0)Tj(c-co)’ 


The  solution  of  (7)  can  be  obtained  by  considering  the 
(1)  equivalent  left-preconditioned  system 


where  Ja  and  c0  are  the  values  of  J  and  c  at  the  minimum, 
respectively,  and  J  is  the  Hessian  matrix  of  the  assimi¬ 
lation  problem. 

Assume  now  for  simplicity  a  linear  framework  and 
minimization  of  the  following  cost  function  with  respect 
to  the  model  state  x  at  the  initial  time  /  =  0  xq  =  c: 


/  =  -(c-c6)tB-1(c-c6) 


I(H„x^d:)TR-'(Hx„-d:) 


(2) 


PTSc  =  PTb  (8) 

with  the  only  requirement  for  the  preconditioner  P  is  to 
have  the  same  range  with  S  (e.g.,  Hayarni  et  al.  2007).  As 
seen  in  (8),  the  best  preconditioner  is  P*  =  SJ-1,  which 
delivers  convergence  in  one  iteration.  In  practice,  how¬ 
ever,  (pseudo)  inversion  of  the  Hessian  is  computa¬ 
tionally  unfeasible,  and  various  approximations  to  P* 
are  considered.  The  major  principle  of  a  preconditioned 
iterative  solver  is  to  implicitly  build  the  solution  of  the 
form 


Here  cb  is  the  first  guess  (background)  model  state  at  1  =  0; 
B  is  its  error  covariance;  summation  is  made  over  obser¬ 
vation  times  tn,  when  the  data  d*  are  available;  R  is  the 
observation  error  covariance;  operators  H„  project  the 
respective  model  states  \„  onto  the  vectors  d*  of  observed 
quantities;  and  the  superscripted  T  stands  for  transposition. 

The  cost  function  (2)  can  be  rewritten  in  terms  of 
deviations  x  =  x  —  Xb  of  the  model  states  from  their 
background  values: 


/ 


x,'B 


l- 


+  X(H„x, 


(3) 


with  H„  =  R  1/2  H  and  d„  =  R  1/2d„  —  H„x*.  Introducing 
the  notation  M"  for  the  model  operator  mapping  c  onto 
x„  and  omitting  tildes  for  convenience  yields  the  explicit 
expression  for  J  in  terms  of  c: 


J  = 


1 

2 


cTB_1c  +  X(H„M”c 

n 


(4) 


Equation  (4)  can  be  represented  in  the  complete 
square  form: 


/  =  !(  Sc  — b)2, 


(5) 


co=c6+Iam(PTS)mPTr0;  r0  =  Scft-b,  (9) 

m 

which  can  be  performed  by  a  variety  of  methods  in¬ 
volving  Krylov  subspaces  tC  (e.g.,  Vuik  et  al.  1996).  In 
particular,  the  well-known  conjugate  gradient  minimi¬ 
zation  algorithm  can  be  viewed  as  a  version  of  the 
Krylov  subspace  technique  with  P  =  S;  that  is,  J  1  is 
replaced  by  the  identity  operator  I. 

The  adjoint -free  variational  methodology  (Yaremchuk 
et  al.  2009)  avoids  the  use  of  the  adjoint  model  M"T  in  the 
construction  of  the  preconditioner  but  employs  explicit 
computation  of  the  columns  of  S  to  invert  J  in  JC  and 
performs  consecutive  orthogonalization  of  the  search 
subspaces  with  respect  to  the  inner  product  induced  by 
the  Hessian  matrix  in  the  control  space  (see  the 
appendix). 

In  principle,  the  availability  of  the  algorithm  for  in¬ 
version  of  J  in  a  sequence  of  J-orthogonal  subspaces  is  a 
necessary  condition  for  the  monotonic  decrease  of  the 
cost  function  with  iterations,  but  it  is  not  sufficient  to 
achieve  the  global  minimum  in  the  linear  case  due  to  a 
possibility  of  the  breakdown  of  the  process.  The  latter 
occurs  when  the  new  search  directions  become  (almost) 
linearly  dependent  on  the  directions  generated  at  the 
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previous  iterations.  In  the  realistic  nonlinear  applications, 
however,  the  number  of  iterations  rarely  exceeds  a  few 
hundred  due  to  computational  constraints.  Therefore,  the 
key  practical  requirement  to  the  algorithm  for  the  itera¬ 
tive  update  of  the  search  directions  is  its  ability  to  quickly 
generate  subspaces  spanned  by  the  leading  projections  of 
the  initial  residual  ro  on  the  eigenvectors  of  the  inverse 
Hessian  matrix,  which  defines  the  solution  copt  =  J  1  STb 
to  (7).  If  the  final  residual  is  sufficiently  small  (i.e.,  the 
model-data  misfit  is  smaller  than  the  observation  error), 
the  solution  is  assumed  to  be  acceptable.  An  ad  hoc 
method  is  to  employ  the  sequence  {M"rm},  as  it  contains 
projections  of  the  residual  rm  at  the  mth  iteration  on  the 
most  persistent  dynamical  modes  of  the  system  and  may 
therefore  capture  a  significant  part  of  copt- 
In  the  present  study  we  employ  the  technique  of 
building  the  search  subspaces  1C  on  the  sequences 
This  approach  provided  acceptable  solutions 
after  100  iterations  and  therefore  proved  to  be  satis¬ 
factory  in  the  experiments  with  WAM. 

The  a4DVAR  algorithm  was  coded  as  the  following 
sequence  of  operations: 

1)  Specify  the  background  initial  conditions  c  and  run 
the  model,  writing  down  model  states  x„  and  their 
deviations  from  the  data  H„x„  —  d*  at  the  observa¬ 
tion  sampling  frequency  (15  min).  Using  the  model 
run  output,  compute  the  cost  function  value  /o  and 
the  auxiliary  vector  Sc. 

2)  Extract  the  m  leading  EOFs  e7,  j  =  1,  m  from  the 
sequence  {x„}  to  form  the  basis  in  the  search  subspace. 

3)  Perturb  the  initial  condition  c,  — >  c  +  ee;-  (e  =  0.001) 
and  perform  (in  parallel)  the  ensemble  of  the  per¬ 
turbed  model  runs,  computing  the  respective  per¬ 
turbed  values  of  Jj  and  Sc,. 

4)  J-orthogonalize  the  search  basis  {e,}  with  respect  to 
at  most  k  bases  obtained  on  the  previous  iterations 
and  compute  the  optimal  correction  5c  to  the  initial 
condition  by  minimizing  J  in  the  subspace  spanned 
by  the  perturbations  (see  the  appendix). 

5)  Update  the  initial  condition  c  <—  c  +  5c  and  compute 
the  updated  values  of  x" ,  J,  and  Sc  by  running  the  model. 

6)  If  the  stopping  criterion  is  satisfied  (either  the 
magnitude  of  the  cost  function  gradient  normalized 
by  its  initial  value  is  less  than  10-3/o  or  the  number  of 
iterations  exceeds  100),  exit.  Otherwise,  go  to  1  to 
start  the  next  iteration. 

b.  Setting  the  experiments 

1)  WAM 


wave  component  with  the  wavenumber  k  =  (kx,  kv)  at 
the  location  x  =  (x,  y)\ 

dl  +  \.(vF)=S(F,x,k,t),  (10) 

at 

where  S  is  the  sum  of  source  functions,  primarily 
composed  of  wind-forced  generation,  dissipation,  and 
redistribution  of  the  wave  spectrum  by  nonlinear  wave- 
wave  interactions  (WAMDI  Group  1988);  V  =  (Vx,  V*} 
stands  for  the  gradient  in  the  horizontal  and  wave- 
number  coordinates;  and  v  is  the  four-component  vector 
of  the  respective  wave-propagation  velocities  depending 
on  the  ambient  current  and  constrained  by  the  disper¬ 
sion  relationship  for  linear  surface  waves: 

cr2  =  gjk|  tanh|kj/7 ,  (11) 

where  cr  is  the  wave  angular  frequency  and  h(x)  is  the 
water  depth. 

Given  the  appropriate  initial/boundary  conditions, 
ambient  current,  and  wind  forcing,  (10)  is  integrated 
numerically  to  produce  the  evolution  of  the  wave  spec¬ 
trum  in  the  space-time  domain  of  the  model. 

An  important  diagnostic  formula  available  in  the 
WAM  package  relates  the  squared  SWH  Q2  with 
the  spectral  density  through  the  following  linear 
relationship: 

Q2(x,t)  =  QN  =  16jJF(x,k,t)dk,  (12) 

k 

where  dk  denotes  the  gridccll  area  in  the  wavenumber 
space  and  summation  is  done  over  the  entire  grid. 

The  model  was  configured  in  the  domain  shown  in 
Fig.  1  with  the  spatial  resolution  of  Sx  =  9  km.  There 
were  mx  =  4412  active  grid  points  in  horizontal  and 
m/c  =  600  grid  points  (24  directions  at  15°  resolution  and 
25  logarithmically  spaced  frequencies  between  0.0314 
and  0.3091  Hz)  in  the  wavenumber  space.  The  total 
length  of  the  state  vector  was  M  =  mx  X  m *  =  2  647  200. 

Distance  between  the  model  states  was  assessed  in 
terms  of  the  correlation  coefficient  C  and  the  normal¬ 
ized  rms  difference  S  between  the  spectra: 


an 


S(F ) 


/  1  2/  ;  N'=N-  (N) 

-|  1/2 

y  k)(f?) 


(13) 


(14) 


WAM  performs  a  time  integration  of  the  balance 
equation  describing  spectral  density  F(x,  k,  t)  for  the 


where  angular  brackets  denote  averaging  in  space,  time, 
and  over  the  wavenumbers.  Similar  coefficients  were 
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calculated  to  assess  the  differences  between  the  scalar 
(SWH)  and  vector  (wind  speed)  fields,  with  averaging 
performed  just  in  space  and  time. 

2)  Cost  function 

The  general  form  of  the  cost  function  used  in  the  data 
assimilation  experiments  was  identical  to  (4)  with  the 
M-dimensional  vector  c  =  N(to)  —  Nb(to )  describing  the 
difference  between  the  gridded  model  state  F(x,  k,  to) 
and  the  background  (first  guess)  state  F),(x,  k,  to)  at  the 
start  of  model  integration  to-  The  first  term  in  (4)  was 
specified  by 

ctB^1c  =  1TX[(I-«2VJ)0c]2  (15) 

X 

and  was  kept  intact  through  all  the  experiments.  In  the 
aforementioned  equation,  W  is  the  regularization 
weight;  I  is  the  mx  X  mx  identity  matrix;  and  the  pa¬ 
rameter  a  defines  the  high-frequency  cutoff  in  horizontal 
coordinates,  which  enforces  spatial  smoothness  in  the 
deviations  of  the  optimal  state  SWH  from  the  back¬ 
ground.  By  setting  a  =  2Sx  throughout  the  experiments, 
SWH  variability  at  spatial  scales  below  a  was  heavily 
penalized.  The  regularization  weight  W  was  chosen  to  be 
inversely  proportional  to  the  squared  mean  of  SWH  in 
the  background  solution  with  the  proportionality  co¬ 
efficient  sx  =  0.01. 

As  seen  in  (15),  the  parameter  a  can  also  be  in¬ 
terpreted  as  the  horizontal  decorrelation  scale  of  the 
SWH  errors:  The  respective  inverse  error  correlation  is 
proportional  to  (I  —  a2V2)2  and  has  the  kernel  that  ex¬ 
ponentially  decays  with  the  distance  between  the  cor¬ 
related  points  at  the  length  scale  of  \\j2h ra  ~  60  km 
(e.g.,  Yaremchuk  and  Smith  2011).  Since  B  1  is  also 
proportional  to  (I  —  a2V2)2,  a  similar  spatial  structure  is 
induced  in  the  spatial  correlations  between  the  spectral 
components  of  c.  Apparently,  our  choice  of  the  decor¬ 
relation  scale  could  be  refined  (see,  e.g,  Waters  et  al. 
2013),  but  we  consider  it  to  be  a  reasonable  approxi¬ 
mation  to  reality  given  the  objectives  of  our  study.  In  the 
spectral  subspace,  (15)  defines  the  inverse  error  co- 
variance  to  have  only  one  linearly  independent  column 
(specified  by  the  components  of  Q ).  As  a  consequence, 
spectral  correlations  at  a  given  point  are  represented  by 
an  m k  X  mu,  correlation  matrix  whose  elements  are  equal 
to  1  (thus  implying  100%  correlation  between  all  the 
spectral  components).  This  assumption  has  been  rou¬ 
tinely  used  in  the  sequential  algorithms  assimilating 
SWH  (e.g.,  Wittmann  and  Cummings  2004). 

The  observation  part  of  the  cost  function  (4)  is  de¬ 
fined  by  the  structure  of  the  observation  operators  H„ 
and  the  respective  error  variances  R„  described  in 
section  2b(4). 


To  compare  the  results  with  the  traditional  OI 
method,  we  used  the  2D  OI  approach  (e.g.,  Wittmann 
and  Cummings  2004;  Waters  et  al.  2013)  in  application 
to  the  SWH  data:  at  the  observation  times,  the  WAM 
state  was  sequentially  updated  by  the  OI  analysis  of  the 
SWH  field,  which  was  projected  onto  the  spectral  com¬ 
ponents  by  multiplying  the  spectrum  at  a  grid  point  by 
the  ratio  of  the  updated  to  predicted  SWH  values.  The 
OI  algorithm  was  configured  with  the  same  background 
error  covariance  B,  R„,  H„  and  using  the  same  true  and 
background  solutions  as  the  a4DVAR  method. 

3)  True  and  background  solutions 

To  perform  the  experiments,  the  true  evolution  of 
the  wave  field  was  generated  by  integrating  WAM  from 
the  state  of  rest  for  10  days  under  Oceanweather  Inc. 
(OWI)  wind  forcing.  The  high-resolution  OWI  winds 
were  developed  by  Oceanweather,  Inc.  using  the 
methodology  of  Cardone  et  al.  (1995, 1996).  The  winds 
were  taken  for  the  period  11-20  September  2011  and 
were  updated  every  hour  during  the  integration.  The 
true  state  Nt(x,  k,  t)  shown  in  Fig.  2  was  picked  from  the 
last  9  h  of  the  model  run  (0000-0900  UTC  20  September 
2011).  Synthetic  data  (described  in  the  next  subsection) 
were  picked  from  the  true  solution  and  then  used  for 
its  reconstruction  through  minimization  of  the  cost 
function  (4). 

Minimization  started  from  the  background  run  of  the 
model,  which  was  obtained  as  follows:  The  true  model 
solution  was  averaged  in  time  and  space  and  the 
resulting  spatially  homogeneous  spectrum  was  used  as 
the  initial  condition  for  the  background  model  run.  The 
run  was  forced  by  the  winds,  which  were  different  from 
those  forcing  the  true  solution.  First,  the  true  winds  were 
horizontally  smoothed  to  mimic  the  errors  typical  for 
reanalysis  winds  from  meteorological  centers  that  are 
usually  available  at  a  coarser  (0.25°-l°)  resolution  and 
have  to  be  interpolated  on  the  fine  resolution  of  grid  of  a 
regional  ocean  model.  In  the  case  considered,  the 
smoothing  was  done  by  the  isotropic  Gaussian  filter  with 
the  half-width  of  25  km.  After  smoothing,  the  winds 
were  rotated  35°  counterclockwise  to  increase  their 
distance  from  the  truth  to  .S’wmci  =  0.67  [(14)].  The  larger 
distance  from  the  true  forcing  was  needed  for  better 
assessment  of  the  observation  impact  on  the  re¬ 
construction  of  initial  conditions,  whose  signature  usu¬ 
ally  persists  for  3-5  h  in  a  typical  wave  model 
integration. 

4)  Synthetic  observations 

In  this  study  two  types  of  data  are  considered:  moored 
observations  of  the  wave  spectra  and  SWH  measure¬ 
ments  from  coastal  HF  radars  and  satellites.  Each  data 
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Fig.  2.  Evolution  of  the  horizontally  averaged  (left)  true  and  (right)  background  spectra. 
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type  had  a  specific  structure  of  the  observation  operator, 
noise,  and  data  generation  procedure. 

Two  tested  mooring  sites  are  shown  in  Fig.  1:  one  was 
positioned  in  the  open  Chukchi  Sea  150  km  northwest  of 
Cape  Hope  and  the  location  of  the  second  mooring 
corresponds  to  the  real  instrument  maintained  in  the 
Bay  of  Kotzebu  from  July  to  December  2007  (Francis 
et  al.  2011).  Observations  d*  by  each  mooring  were 
generated  by  multiplying  the  true  spectrum  at  any  mo¬ 
ment  by  the  random  factor  1  +  sv,  where  v  is  the  white 
noise  with  unit  variance  and  e,„  =  0.01.  The  observa¬ 
tional  error  covariance  matrices  R  for  both  moorings 
were  diagonal  with  time-independent  diagonal  elements 
equal  to  (sN,)2.  The  respective  observational  operators 
H  were  time  independent  and  picked  the  time-varying 


WAM  spectra  every  15  min  from  the  grid  point  nearest 
to  the  buoy  location,  providing  4m*  =  2400  observations 
per  hour. 

It  should  be  noted  that  directional  buoys  provide 
only  a  few  moments  of  the  full  spectrum  and  therefore 
constrain  much  smaller  degrees  of  freedom  of  a  model 
than  m*  (usually,  2-3  times  the  number  of  frequencies, 
i.e.,  50-80  in  our  case).  These  moments,  however,  are 
linear  functions  in  the  model  state,  and  their  calculation 
in  terms  of  spectral  components  does  not  principally 
differ  from  Eq.  12.  Assimilation  of  buoy  moments 
would,  of  course,  be  a  better  approach  in  practice  [like 
extending  the  method  of  Voorrips  et  al.  (1997)  to  all  the 
measured  components  of  the  directional  spectrum],  but 
in  context  of  our  study,  this  would  fall  somewhere  in 
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between  assimilating  the  full  spectrum  and  its  first  mo¬ 
ment  (SWH)  and  would  not  be  a  useful  addition  to  our 
analysis.  We  should  also  keep  in  mind  that  a  user  can 
assimilate  postprocessed  buoy  data  in  the  form  of  the 
full  spectrum  computed  by  the  estimator  available  in  the 
software  provided  with  an  instrument. 

SWH  observations  were  simulated  by  integrating  the 
true  spectrum  held  (Eq.  12)  into  the  apexes  of  the  grid 
cell  containing  a  high-frequency  radar  (HFR)  observa¬ 
tion  point  followed  by  linear  interpolation  onto  that 
point.  After  that,  the  SWH  value  was  contaminated  by 
random  noise  with  the  rms  variance  of  30  cm.  HF  radar 
observation  points  were  located  along  the  beams  of  the 
radar  considered  for  placement  at  the  Cape  Hope  site  to 
monitor  wave  conditions  north  of  the  Bering  Strait 
(Fig.  1).  The  above-described  HFR  observation  opera¬ 
tor  computed  SWH  values  along  the  25  beams  every 
15  min.  providing  information  to  535  model  grid  points 
within  the  sector  shown  in  Fig.  1  (2140  observations 
per  hour). 

Synthetic  satellite  observations  of  sea  surface  rough¬ 
ness  provided  SWH  data  along  the  Environmental  Sat¬ 
ellite  ( Envisat )  tracks  shown  in  Fig.  1  with  9-km 
discretization  (55  and  73  points  for  tracks  A  and  B,  re¬ 
spectively).  These  data  were  assumed  instantaneous  and 
satellite  passage  occurred  for  both  tracks  after  2h  of 
model  integration.  The  respective  observation  operator 
was  similar  to  the  one  used  for  HFR,  except  that  it 
picked  SWH  values  at  the  sequence  of  WAM  grid  points 
closest  to  the  sampling  points  along  the  tracks  (i.e.,  no 
spatial  interpolation  was  used).  Satellite  SWH  obser¬ 
vations  were  contaminated  similarly  to  HFRs  with  the 
rms  error  variance  of  30  cm. 

3.  Results 

Synthetic  observations  of  SWH  and  wave  spectra 
were  assimilated  into  WAM  using  the  adjoint-free  var¬ 
iational  technique.  The  model  was  constrained  by  data 
during  the  first  3h  of  model  integration  and  then  in¬ 
tegrated  for  6  h  to  assess  the  improvement  of  the  fore¬ 
cast  skill.  Performance  of  the  method  was  quantified  by 
calculating  C  (Eq.  13)  and  S  (Eq.  14)  between  the  op¬ 
timized  and  true  solutions.  These  quantities  were  com¬ 
puted  with  time  averaging  over  three  time  intervals:  0-3  h 
(assimilation  period),  and  two  forecast  periods  of  3-6  h 
and  6-9  h.  Five  ensemble  members  were  used  for  span¬ 
ning  the  search  subspace  on  every  iteration.  More  de¬ 
tails  on  building  and  orthogonalization  of  the  ensemble 
members  are  given  in  the  appendix.  The  3-h  ensemble 
model  runs  were  executed  in  parallel  and  required  62  s 
of  wall  time  per  a4DVAR  iteration  on  five  processors 
of  the  2.3-GHz  cluster  (1-1.5  h  per  experiment).  The 


iterations  iterations 


Fig.  3.  Convergence  of  the  a4DVAR  algorithm  in  the  assimila¬ 
tion  experiments  with  M1B  and  HF  data,  (left)  Cost  function  and 
(right)  its  gradient  are  normalized  by  their  values  at  the  start  of 
assimilation. 

sequential  OI  algorithm  assimilating  HFR  data  was  ex¬ 
ecuted  in  74  s  on  a  single  processor. 

A  series  of  OI  and  a4DVAR  experiments  were  con¬ 
ducted,  involving  assimilation  of  the  data  from  five 
sources  and  their  combinations:  high-frequency  radar  at 
Point  Hope  (denoted  by  HF),  two  moorings  (a4DVAR 
analyses  only,  locations  shown  in  Fig.  1),  and  two  En¬ 
visat  tracks  (A  and  B;  Fig.  1).  For  comparison  purposes, 
we  conducted  similar  experiments  with  the  OI  method 
assimilating  only  SWH  data  from  satellites  and/or  HF 
radar,  which  are  abbreviated  oHFA(B)  and  oHF,  re¬ 
spectively,  in  the  description  below.  With  the  exception 
of  satellite  tracks,  all  a4DVAR  assimilation  experiments 
demonstrated  significant  improvement  of  the  model 
state  in  terms  of  its  proximity  to  the  true  solution.  The 
stopping  criterion  for  optimization  was  the  reduction  of 
the  cost  function  gradient  1000  times,  which  usually 
occurred  after  80-100  iterations.  By  that  time  the  value 
of  J  was  typically  reduced  2-3  times  (Fig.  3).  Maps  of 
deviations  from  the  truth  of  the  spatially  averaged 
background  and  optimized  spectra  at  t  =  0  are  shown  in 
Fig.  4.  In  most  of  the  a4DVAR  experiments,  the  initial 
error  has  been  reduced  to  the  values  compatible  with  the 
wind  forcing  errors.  The  only  exceptions  were  the  results 
of  optimal  interpolation  (Fig.  4b)  and  of  the  a4DVAR 
assimilation  of  SWH  data  from  a  single-satellite  track 
(not  shown):  in  these  cases  the  optimized  spectrum  was 
only  slightly  different  from  the  one  produced  by  the 
background  solution.  For  the  OI  case  such  a  small 
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Fig.  4.  Absolute  difference  between  the  horizontally  averaged  true  spectrum  at  t  =  0  and 
(a)  background  spectrum,  (b)  oHFA-optimized,  (c)  MlA-optimized,  and  (d)  HFA-optimized 
spectra. 


correction  can  be  explained  by  the  fact  that  SWH  data 
are  weakly  constrained  by  dynamics  and  can  barely  af¬ 
fect  the  shape  and  location  of  the  spectra  because  they 
provide  information  only  on  their  mean  magnitude  at  a 
geographical  position.  Small  spectral  improvement  of 
the  a4DVAR  experiments  with  a  single-satellite  track 
could  be  attributed  to  the  small  amount  of  data  (55  SWH 
observations).  As  a  consequence,  the  cost  function  is 
dominated  by  the  regularization  term,  which  implies 
100%  correlations  in  spectral  space  and  is  therefore 
capable  of  adjusting  only  the  spectral  magnitude. 

These  properties  of  the  abovementioned  assimilated 
solutions  translate  into  their  lower  spectral  forecast  skill 
shown  in  Table  1,  which  also  includes  spectral  errors  from 
the  other  assimilation  experiments.  Abbreviations  in  the 
left  column  of  Table  1  correspond  to  the  types  of  data 
used  in  the  experiment  (e.g.,  HFA  corresponds  to  as¬ 
similation  of  the  HF  data  and  SWH  data  from  the  Envisat 
track  A).  As  seen,  6420  observations  by  the  HF  radar 
provide  the  largest  impact  on  the  improvement  of  the 
wave  field  during  the  assimilation  period:  the  correlation 
coefficient  C03  increases  from  0.47  to  0.77,  whereas  S® 
drops  from  0.89  to  0.65.  The  OI  method  provides  a  rela¬ 
tively  low  skill  improvement  similar  to  assimilation  of  a 
single-satellite  track  (cf.  row  1  and  rows  5,  13,  14). 


Direct  measurement  of  the  spectra  by  a  single  moor¬ 
ing  (7200  observation  points,  rows  6,  9)  also  provide 
only  a  moderate  increase  of  C®  to  0.52  and  a  decrease  of 
S®  to  0.86.  This  can  be  partly  explained  by  the  fact  that 
assimilated  spectra  occupy  a  small  part  of  spectral  do¬ 
main  (at  most  15%-20%;  Fig.  2).  As  a  consequence,  the 
effective  number  of  observations  with  useful  (nonzero) 
information  on  the  state  of  the  wave  field  should  be 
reduced  5-7  times  down  to  —1500  data  points  on  the 
total,  which  is  compatible,  by  the  way,  to  assimilating  3- 
7  spectral  moments.  Besides,  mooring  data  do  not  pro¬ 
vide  any  information  on  the  spatial  variability  of  the 
spectra,  which  appears  to  be  crucial  for  the  successful 
recovery  of  the  true  state. 

In  that  respect,  it  is  remarkable  that  adding  much  less 
numerous  satellite  data  to  moored  spectral  observations 
improves  the  performance  of  the  assimilation  system 
considerably.  Combining  moored  and  satellite  data 
provides  30%-40%  growth  of  the  correlation  co¬ 
efficients  and  a  20%-25%  drop  of  the  normalized 
standard  deviations  from  the  true  spectrum  (cf.  rows  6 
and  9  with  rows  7-8  and  10-11,  respectively).  At  the 
same  time,  satellite  SWH  data  do  not  add  much  new 
information  to  that  containing  in  HFR  observations  (cf. 
rows  2  and  3-4),  which  monitor  the  same  integral 
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Table  1.  Normalized  S  and  C  between  the  optimized  and  true 
solutions  for  the  experiments  with  various  types  of  data.  Subscripts 
03, 36,  and  69  correspond  to  time-averaging  ranges  of  0-3,  3-6,  and 
6-9  h  of  model  integration. 


No. 

Expt 

Co3 

u>03 

C36 

S36 

Cb9 

^69 

1 

BG 

0.47 

0.89 

0.48 

0.87 

0.59 

0.80 

2 

HF 

0.77 

0.65 

0.72 

0.69 

0.71 

0.70 

3 

HFA 

0.75 

0.66 

0.72 

0.68 

0.70 

0.70 

4 

HFB 

0.76 

0.64 

0.76 

0.65 

0.75 

0.67 

5 

oHFA 

0.48 

0.87 

0.49 

0.86 

0.59 

0.79 

6 

Ml 

0.53 

0.85 

0.47 

0.87 

0.50 

0.86 

7 

MIA 

0.70 

0.71 

0.70 

0.71 

0.68 

0.73 

8 

M1B 

0.71 

0.70 

0.70 

0.71 

0.69 

0.72 

9 

M2 

0.52 

0.87 

0.50 

0.87 

0.57 

0.86 

10 

M2A 

0.65 

0.76 

0.69 

0.72 

0.67 

0.73 

11 

M2B 

0.69 

0.75 

0.69 

0.72 

0.68 

0.73 

12 

M1M2 

0.71 

0.73 

0.67 

0.75 

0.68 

0.73 

13 

A 

0.47 

0.89 

0.48 

0.88 

0.57 

0.79 

14 

B 

0.48 

0.88 

0.48 

0.87 

0.58 

0.78 

quantity  for  the  whole  assimilation  period  (3h)  and 
cover  a  significant  part  of  the  model  domain  (Fig.  1). 

The  importance  of  the  spatial  coverage  by  observa¬ 
tions  is  confirmed  by  the  result  of  the  experiment  with 
assimilation  of  the  spectra  from  two  moorings:  The 
values  of  C  and  S  in  this  case  demonstrate  a  considerable 
improvement  and  become  compatible  (row  12  in  Table 
1)  with  those  achieved  with  the  joint  assimilation  of 
spectra  from  the  single  mooring  and  satellite  SWH  data 
(rows  7-8  and  10-11). 

Normalized  rms  deviation  term  S  from  the  true  solu¬ 
tion  as  a  function  of  time  is  shown  in  Fig.  5.  It  is  note¬ 
worthy  that  the  moderate  reduction  of  S  acquired  after 
assimilating  Mi  data  is  completely  lost  after  2h  of  in¬ 
tegrating  the  model  without  data  constraints:  after  /  =  5 
hours,  deviation  from  the  truth  becomes  even  larger 
than  that  of  the  OI  solution,  which  is  very  close  to  the 
respective  deviation  of  the  background.  The  situation  is 
much  improved  with  satellite  data  being  taken  into  ac¬ 
count  (dashed  lines  in  Fig.  5):  the  value  of  S  significantly 
drops  and  becomes  compatible  with  the  values  achieved 
by  assimilating  HFR  observations.  Complementing 
HFR  data  with  satellite  observations  does  not  affect  the 
mean  value  of  S  during  the  assimilation  period,  but  it 
has  a  noticeable  impact  on  the  forecast  skill  of  the  model 
(quantified  by  S )  up  to  the  end  of  integration.  Apparent 
convergence  of  all  the  curves  by  hour  9  is  the  conse¬ 
quence  of  the  wind  forcing,  which  tends  to  attract  model 
trajectories  to  a  state  whose  error  is  compatible  with  that 
of  the  wind  (SWmd  =  0.67). 

The  time  dependence  of  the  spatially  averaged  cor¬ 
relation  coefficients  C(t)  between  the  true  and  optimized 
spectra  (Fig.  6)  demonstrates  similar  behavior:  assimi¬ 
lation  of  the  M  |  spectra  provides  a  slight  improvement  of 


t,  hours 

Fig.  5.  Normalized  rms  deviations  from  the  truth  of  the  oHFA, 
HF-optimized,  and  Ml-optimized  spectra.  Dashed  lines  corre¬ 
spond  to  augmenting  the  a4DVAR-assimilated  data  with  satellite 
SWH  observations  from  satellite  passing  through  track  B  (experi¬ 
ments  HFA  and  HFB).  Thick  gray  line  corresponds  to  the  M1-M2- 
optimized  a4DVAR  solution.  Vertical  line  marks  the  end  of 
assimilation  period. 

the  correlation  (cf.  gray  and  thick  black  curve  in  Fig.  6), 
which  persists  for  2h  of  model  integration.  Adding  sat¬ 
ellite  SWH  data  from  track  B  instantly  provides  in¬ 
formation  on  the  higher  spectral  densities  in  the 
northern  part  of  the  Chukchi  Sea  (Fig.  1)  and  rises  the 
correlation  to  0.65-0.7  (dashed  gray  line  in  Fig.  6). 

Inspection  of  Table  1  also  shows  that  information 
from  track  A  also  increases  the  efficiency  of  assimilating 
spectra  from  moorings  but  to  somewhat  lesser  extent. 
This  phenomenon  can  be  partly  explained  by  the  fact 
that  track  A  does  not  cover  the  region  of  the  highest 
SWH  and  therefore  provides  less  information  on  the 
magnitude  of  spatial  variability  of  the  wave  field.  Simi¬ 
larly,  assimilation  of  the  M2  data  appears  to  be  slightly 
less  efficient  than  M1;  which  can  be  partly  attributed  to 
the  M2  position  at  the  periphery  of  the  domain. 

Table  1  also  indicates  that  on  a  regional  scale,  in¬ 
stantaneous  Envisat  observations  cannot  provide  a  sig¬ 
nificant  improvement  to  the  background  state  if  they  are 
not  accompanied  by  in  situ  measurements.  At  the  same 
time,  satellite  data  become  quite  valuable  in  com¬ 
plementing  observations  if  the  wave  conditions  are 
measured  by  a  single  mooring. 

The  forecast  errors  provided  by  the  HFA  data  as¬ 
similation  using  OI  and  a4DVAR  techniques  and 
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t,  hours 

FIG.  6.  As  in  Fig.  5,  but  for  correlations. 

averaged  over  the  period  of  3-9  h  are  compared  in  Fig.  7 
in  terms  of  the  horizontal  distributions  of  the  SWH, 
peak  period,  and  wave  direction  errors,  that  The 
a4DVAR  technique  provides  30%-50%  better  forecast 
skill  in  terms  of  the  SWH  (0.28  vs  0.37  m)  and  peak 
period  (0.63  vs  0.91  s).  Although  discrepancies  in  the 
peak  period  near  the  southern  and  eastern  boundaries 
are  comparable  in  both  solutions,  the  a4DVAR  method 
demonstrates  a  significant  advantage  over  OI  in  the 
northern  Chukchi  Sea  and  south  of  Cape  Hope,  result¬ 
ing  in  approximately  10-cm  smaller  SWH  errors 
throughout  the  entire  domain.  A  local  maximum  in  the 
a4DVAR  peak  period  errors  is  also  observed  southwest 
of  Cape  Hope  (Fig.  7b)  that  can  be  partly  explained  by  a 
sharper  gradient  in  the  peak  period  field  of  the  true  so¬ 
lution  (Fig.  7a). 

The  OI  solution  demonstrates  a  slightly  better  skill  in 
forecasting  the  wave  direction  (the  mean  difference  of 
13.1°  vs  16.9°).  However,  in  the  OI  assimilation  exper¬ 
iments  with  other  types  of  SWH  data  this  number 
varied  within  13°-13.3°  and  was  quite  close  to  the  re¬ 
spective  characteristic  (13.2°)  of  the  background 
solution. 

In  general,  our  experiments  have  shown  that  the  OI 
method  tends  to  improve  the  amplitude  of  the  spectrum 
and  that  it  has  only  a  slight  impact  on  its  shape  and 
position  in  the  frequency-direction  coordinates.  In 
contrast,  the  a4DVAR  technique  is  capable  of  improv¬ 
ing  these  characteristics  as  well,  since  it  performs  opti¬ 
mization  along  the  most  persistent  dynamical  modes  of 


180  176  172  168  164  160W 


Fig.  7.  Time-averaged  (3-9  h)  values  of  SWH  (contours,  m), 
peak  period  (shading,  s),  and  wave  direction  (arrows)  of  (a)  the 
true  state  and  the  respective  absolute  differences  of  the  (b)  HF- 
optimized  and  (c)  oHF  solutions  from  the  true  state.  Domain- 
averaged  values  of  the  fields  are  shown  on  the  right. 

the  governing  (10).  This  important  property  of  the 
a4DVAR  algorithm  provides  a  significantly  better  ap¬ 
proximation  of  the  true  solution  and  an  improved 
forecast  skill. 
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4.  Conclusions  and  discussion 

In  this  study,  the  adjoint-free  variational  data  as¬ 
similation  algorithm  of  Yaremchuk  et  al.  (2009)  has 
been  applied  to  the  fully  nonlinear  regional  WAM 
configured  in  the  Chukchi  Sea.  Performance  of  the 
method  was  tested  in  a  series  of  twin-data  experiments 
with  an  objective  to  assess  the  efficiency  of  various 
wave-observing  platforms  in  reconstruction  of  the 
surface  waves  on  a  regional  scale.  Results  indicate  that 
an  HF  radar  monitoring  of  the  SWH  field  is  more  ef¬ 
ficient  than  a  mooring  measuring  the  wave  spectrum 
at  a  single  point.  However,  augmenting  spectral  ob¬ 
servations  from  the  mooring  with  satellite  SWH  data 
boosts  the  system’s  efficiency  to  the  level  of  an  HF 
radar,  which  does  not  gain  as  much  new  information 
from  satellite  data. 

Application  of  the  adjoint-free  variational  method  to 
optimization  of  the  initial  conditions  is  a  new  de¬ 
velopment  in  oceanic  wave  modeling  that  provides  a 
pathway  to  bypass  the  burden  of  development  and 
maintenance  of  the  tangent  linear  and  adjoint  code. 
Besides,  the  adjoint-free  approach  is  insusceptible  to 
intrinsic  instabilities  arising  in  the  adjoint  model  when 
the  parent  model  solution  is  contains  nonlinear  in¬ 
stabilities.  The  a4DVAR  method  is  based  on  probing 
the  shape  of  the  cost  function  by  a  sequence  of  ensem¬ 
bles  updated  in  the  iterative  process.  In  this  study  we 
utilized  the  simplest  possible  technique  of  generating 
the  ensemble  members  by  the  EOF  analysis  of  model 
solutions.  Although  our  experiments  demonstrated  a 
reasonably  good  performance  of  the  method,  more  ad¬ 
vanced  schemes,  such  as  building  the  search  subspaces 
(spanned  by  the  ensemble  members)  on  the  cost  func¬ 
tion  gradients,  could  be  developed  to  achieve  faster 
convergence. 

An  extremely  important  issue  of  the  a4DVAR  de¬ 
velopment  in  application  to  wave  modeling  is  optimi¬ 
zation  of  the  wind  forcing,  whose  uncertainties  provide 
the  major  contribution  to  model  errors  at  time  scales 
larger  than  a  few  hours.  One  possible  solution  is  to 
employ  wind  error  statistics  and  augment  the  search 
subspaces  (ensembles)  by  new  dimensions  spanned  by 
the  leading  EOFs  of  the  wind  error  fields.  However, 
extension  of  the  a4DVAR  technique  to  wind  correction 
capability  is  a  separate  theoretical  issue  lying  beyond  the 
scope  of  the  present  paper. 

An  interesting  finding  of  the  experiments  with  syn¬ 
thetic  data  is  the  important  role  of  satellite  data  (spo¬ 
radic  in  the  setting  considered)  in  complementing 
spectral  observations  at  a  single  or  sparsely  located 
moorings.  At  the  same  time,  the  revealed  inadequacy 
of  assimilating  solely  Envisat  data  (rows  13  and  14  in 


Table  1)  is  apparently  an  artifact  of  the  space-time  size 
of  the  domain  and  substantial  wind  forcing  errors.  In  a 
series  of  experiments  with  .Swmd  =  0.1,  satellite  data  pro¬ 
vided  much  better  improvement  of  the  initial  conditions, 
although  not  at  the  levels  of  HFR  observations. 

In  the  experiments  we  did  not  consider  all  types  of 
data  available  for  regional  modeling.  In  particular,  di¬ 
rectional  spectra  (with  a  coarser  resolution  and  larger 
uncertainties)  can  be  obtained  from  high-frequency  ra¬ 
dars  (Graber  and  Heron  1997;  Hisaki  2005)  and  space- 
borne  (Tison  et  al.  2008;  Ren  et  al.  2010)  platforms. 
Simplifying  assumptions  have  been  also  made  in  con¬ 
figuring  WAM:  the  background  currents  (—0.2- 
0.4 ms-1)  were  assumed  to  be  negligible  in  comparison 
with  the  wave  velocities  (— 5-10 ms-1)  of  the  major 
spectral  peak  and  stationary  conditions  at  the  open 
boundaries  were  specified  due  to  relatively  short  period 
of  model  integration.  These  assumptions  can  be  relaxed 
in  the  experiments  with  real  data  spanning  longer  time 
periods  in  smaller  domains. 

The  4D  variational  (either  adjoint  or  adjoint  free) 
methods  of  data  assimilation  require  multiple  model 
runs  and  are  certainly  more  expensive  than  sequential 
methods  (e.g.,  based  on  optimal  interpolation),  which 
update  the  model  state  in  the  process  of  single  in¬ 
tegration.  The  4DVAR  methods  do,  however,  deliver 
better  forecast  skill  and  are  therefore  competitive  with 
sequential  algorithms  when  the  computational  cost  is 
not  a  primary  issue.  In  our  case,  one  assimilation  run  was 
performed  at  the  expense  of  approximately  500  direct 
model  runs,  which  is  compatible  to  the  cost  of  a  typical 
adjoint  4DVAR  run.  The  major  advantage  of  the  pro¬ 
posed  adjoint-free  technique  is  the  absence  of  necessity 
to  develop  and  maintain  the  adjoint  and  tangent  linear 
codes  for  the  assimilative  model. 

Results  of  the  present  study  can  be  used  in  the  design  of 
the  wave  monitoring  systems  supported  by  the  variational 
data  assimilation  into  state-of-the-art  wave  models  (Booij 
et  al.  1996;  Tolman  and  the  WAVEWATCH  III 
Development  Group  2014)  lacking  the  comprehensive 
adjoint  and  tangent  linear  codes.  Nowadays,  de¬ 
velopment  of  such  systems  becomes  particularly  im¬ 
portant  in  the  Arctic  seas,  which  experience  a  significant 
surge  in  SWH  variability  associated  with  ongoing  ice 
retreat  induced  by  climate  change  (e.g.,  Francis  and 
Atkinson  2012;  Thomson  and  Rogers  2014). 
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APPENDIX 

Hessian  Projection  and  Orthogonalization 

In  the  present  study,  each  iteration  of  the  a4DVAR 
minimization  starts  with  the  reference  model  run  from 
the  initial  condition  c  obtained  at  the  previous  iteration. 
This  run  produces  the  a  sequence  of  model  states  subject 
to  principal  orthogonal  decomposition  analysis.  The 
leading  modes  e„,  n  =  1,  . . . ,  N  (N  =  5)  of  the  analysis 
are  utilized  to  perturb  the  reference  initial  condition 
c  — >  c  4-  ee„,  to  execute  N  parallel  model  runs,  and  to 
compute  the  vectors  z„  =  Se„,  whose  mutual  scalar 
products  z^Z(t  =  e'^  Jc/f  can  be  used  to  perform  Gram- 
Schmidt  orthogonalization  in  the  control  space  (span¬ 
ned  by  ek)  with  respect  to  the  inner  product  induced  by 
the  Hessian  matrix  J. 

The  perturbation  vectors  5c„  =  ee„  span  the  search 
subspace,  which  can  be  parameterized  by  the  N- 
dimensional  vector  a  such  that 

5c  =  E  a,  (Al) 

where  e„  constitute  the  columns  of  the  N  X  M  matrix  E 
projecting  the  Hessian  on  the  search  subspace.  Multi¬ 
plying  (7)  by  Et,  and  taking  (Al)  into  account,  yields  the 
projection  of  the  normal  equations  on  search  subspace: 

ETSTSEa  =  ETSTb.  (A2) 

To  solve  (A2)  it  is  necessary  to  compute  both  the  Hes¬ 
sian  projection  ETSTSE  and  the  result  of  action  of  its 
square  root  ETST  on  b.  These  computations  are 
straightforward. 

The  matrix  SE  is  computed  explicitly  columnwise  by 
taking  the  differences  between  the  perturbed  and  un¬ 
perturbed  residuals  S(c  +  ee„)  —  Sc  =  eSe„.  Dividing  the 
resulting  N  vectors  z„  =  eSe„  by  e,  the  columns  of  SE  are 
obtained.  Explicit  computation  of  ETSTSE  =  (SE)TSE 
is  straightforward:  one  has  to  compute  N(N  +  l)/2 
inner  products  of  all  the  possible  pairs  of  z„  (e.g., 
Zupanski  2005). 

In  the  similar  manner,  the  perturbation  vectors  cn  can 
be  orthogonalized  with  respect  to  the  inner  product  in¬ 
duced  by  J  if  the  vectors  z„  obtained  on  the  previous 
iterations  are  kept  in  memory. 

To  assess  ETSTb,  consider  a  small  perturbation  of  /, 
induced  by  e„: 


8Jn  =  ^(5c»sTs5c„  ~  ScJSTb  -  brS5c„  +  bTb)  -  ^ 
=  ^5cTSTS5cn-5cTSTb. 

(A3) 

The  first  term  on  the  rhs  of  (A3)  is  negligible  because 
5c„  =  een  and  e  ^  1.  As  a  consequence, 

8Jn  =  sej  STb.  (A4) 

Collecting  8Jn  into  a  single  vector  5J  and  dividing  by 
e  yields  the  rhs  of  (A2): 

-5J/e  =  ETSTb.  (A5) 

As  a  result,  the  system  of  normal  equations  in  the  search 
subspace  can  be  solved  at  a  relatively  low  cost  if  the 
vectors  z„  and  the  respective  perturbations  of  the  cost 
function  8Jn  are  available.  In  the  a4DVAR  algorithm, 
these  quantities  are  computed  independently  during  the 
parallel  runs  of  the  perturbed  model. 
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